/* Copyright 2012 Tobias Marschall
 *
 * This file is part of CLEVER.
 *
 * CLEVER is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * CLEVER is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with CLEVER.  If not, see <http://www.gnu.org/licenses/>.
 */

#ifndef SPLITALIGNER_H_
#define SPLITALIGNER_H_

#include <iostream>
#include <cstddef>
#include <vector>
#include <memory>
#include <limits>
#include <iomanip>

#include <bamtools/api/BamAux.h>

class SplitAligner {
public:
	typedef enum { LEFT_ANCHOR, RIGHT_ANCHOR } anchor_type_t;
	typedef enum { NO_ALIGNMENT, REGULAR_ALIGNMENT, SPLIT_INSERTION, SPLIT_DELETION } split_alignment_type_t;
	typedef enum { NONE, LEFT, ABOVE, ABOVE_LEFT } backtrace_t;
	typedef enum { LINEAR, AFFINE } gap_costs_t;
	typedef int cost_t;
	/** Type for information on how far a match could be extended. */
	typedef struct anchor_alignment_t {
		size_t total_query_length;
		/** Number of columns that have been computed. */
		size_t column_count;
		/** Number of rows that have been computed. */
		size_t row_count;
		/** Start position of alignment on the reference. */
		size_t ref_start_pos;
		/** Is this a left or a right anchor alignment? */
		anchor_type_t type;
		/** A vector of 1 + 2*max_indels diagonals. The cell (i,j) of a "regular" DP matrix
		 *  is located at diagonals[i-j+max_indels][j]. */
		std::vector<std::vector<cost_t> > diagonals;
		/** The minimum cost for each column. */
		std::vector<cost_t> column_minima;
		/** The minimum cost for each row. */
		std::vector<cost_t> row_minima;
		/** Backtrace table to reconstruct alignment. */
		std::vector<std::vector<backtrace_t> > backtrace;
		/** Indel length table used for affine gap costs to quickly determine length of indel starting from respective cell. */
		std::vector<std::vector<int> > indel_length_table;
		anchor_alignment_t(size_t total_query_length, size_t ref_start_pos, size_t max_indels, anchor_type_t type, cost_t cost_threshold, bool init_indel_length_table) : 
			total_query_length(total_query_length),
			column_count(0),
			row_count(0),
			ref_start_pos(ref_start_pos),
			type(type),
			diagonals(1+2*max_indels, std::vector<cost_t>(total_query_length+1,cost_threshold+1)),
			column_minima(total_query_length+1,cost_threshold+1),
			row_minima(total_query_length+max_indels+1,cost_threshold+1),
			backtrace(1+2*max_indels, std::vector<backtrace_t>(total_query_length+1,NONE)),
			indel_length_table((init_indel_length_table?(1+2*max_indels):0), std::vector<int>(total_query_length+1,-1)) {}
	} anchor_alignment_t;

	typedef struct split_alignment_t {
		split_alignment_type_t type;
		int length;
		// positions right after the breakpoints in reference ...
		int reference_breakpoint;
		// ... and in query
		int query_breakpoint;
		int phred_score;
		int mismatch_phred_score;
		std::vector<BamTools::CigarOp> cigar;
		split_alignment_t() : type(NO_ALIGNMENT), length(-1), reference_breakpoint(-1), query_breakpoint(-1), phred_score(-1), mismatch_phred_score(-1) {}
		split_alignment_t(split_alignment_type_t type, int length, int reference_breakpoint, int query_breakpoint, int phred_score, int mismatch_phred_score) : type(type), length(length), reference_breakpoint(reference_breakpoint), query_breakpoint(query_breakpoint), phred_score(phred_score), mismatch_phred_score(mismatch_phred_score) {}
	} split_alignment_t;

	typedef struct interchrom_split_alignment_t {
		// positions right after the breakpoints in reference ...
		int reference_breakpoint1;
		int reference_breakpoint2;
		int phred_score;
		int mismatch_phred_score;
		// phred-scaled cost incurred by the split
		int split_cost;
		std::vector<BamTools::CigarOp> cigar1;
		std::vector<BamTools::CigarOp> cigar2;
		interchrom_split_alignment_t(int reference_breakpoint1, int reference_breakpoint2, int phred_score, int mismatch_phred_score, int split_cost) : reference_breakpoint1(reference_breakpoint1), reference_breakpoint2(reference_breakpoint2), phred_score(phred_score), mismatch_phred_score(mismatch_phred_score), split_cost(split_cost) {}
	} interchrom_split_alignment_t;

	typedef struct whole_alignment_t {
		int start_position;
		int end_position;
		int phred_score;
		int mismatch_phred_score;
		std::vector<BamTools::CigarOp> cigar;
		whole_alignment_t() : start_position(-1), end_position(-1), phred_score(-1), mismatch_phred_score(-1) {}
	} whole_alignment_t;

private:
	// type of gap costs (linear vs. affine)
	gap_costs_t gap_costs;
	// offset to be subtracted from ASCII values in PHRED strings
	int phred_offset;
	// "PHRED like" penalty for an indel, if gap_costs==AFFINE, then this value contains the gap open costs
	int indel_cost;
	// "PHRED like" cost for extending a gap (only used when gap_costs==AFFINE)
	int gap_extend_cost;
	// "PHRED like" penalty for splitting a read
	int split_cost;
	// "PHRED like" penalty for splitting a read interchromosomally (if that's a word ;-)
	int interchrom_split_cost;
	// "PHRED like" threshold for cost sum, when threshold is exceeded, no
	// alignment is reported
	int cost_threshold;
	// The maximal number of allowed indels (and thus the width of the band for
	// the banded alignment).
	int max_indels;
	// Report all split alignments with a PHRED score difference of at most this
	// value to the best split alignment.
	int secondary_aln_phred_diff;
	// minimal length of anchor on both sides of a split alignment
	int min_anchor_length;
	// if true, only M characters (instead of X and =) are used for cigar strings
	bool simple_cigars;
	// report rightmost alignments (instead of leftmost)
	bool rightmost;
	bool allow_soft_clipping;
	// cost for soft clipping a read
	int soft_clip_open_cost;
	// cost for "extending" a soft clip; i.e., softclipping k characters from a read 
	// will cost soft_clip_open_cost+k*soft_clip_extend_cost
	int soft_clip_extend_cost;

	typedef struct split_alignment_comparator_t {
		bool rightmost;
		bool operator()(const split_alignment_t& a, const split_alignment_t& b) { 
			if (a.phred_score != b.phred_score) return a.phred_score < b.phred_score;
			if (rightmost) {
				return a.reference_breakpoint > b.reference_breakpoint;
			} else {
				return a.reference_breakpoint < b.reference_breakpoint;
			}
		}
		split_alignment_comparator_t(bool rightmost) : rightmost(rightmost) {}
	} split_alignment_comparator_t;

	typedef struct interchrom_split_alignment_comparator_t {
		bool rightmost;
		bool operator()(const interchrom_split_alignment_t& a, const interchrom_split_alignment_t& b) { 
			if (a.phred_score != b.phred_score) return a.phred_score < b.phred_score;
			return a.reference_breakpoint1 < b.reference_breakpoint1;
		}
		interchrom_split_alignment_comparator_t() {}
	} interchrom_split_alignment_comparator_t;

	template <typename IteratorType>
	void append_cigar(const IteratorType& begin, const IteratorType& end, std::vector<BamTools::CigarOp>* target) const {
		for (IteratorType it = begin; it != end; ++it) {
			if ((target->size() > 0) && (target->at(target->size()-1).Type == it->Type)) {
				target->at(target->size()-1).Length += it->Length;
			} else {
				target->push_back(*it);
			}
		}
	}

	/** Compute the reverse of the cigar string of the alignment starting at given diagonal and column,
	 *  and append it to the given target vector. Adds costs incurred by mismatches to *mismatch_phred_score. */
	void compute_reverse_cigar(const anchor_alignment_t& anchor, int diag, int col, std::vector<BamTools::CigarOp>* target, int* mismatch_phred_score) const;
	void anchor_to_cigar(const anchor_alignment_t& anchor, int diag, int col, std::vector<BamTools::CigarOp>* cigar, int* mismatch_phred_score) const;
	/** Prepare an deletion-like instance of split_alignment_t and push it to a vector. */
	void push_deletion(const anchor_alignment_t& left, const anchor_alignment_t& right,  int left_ref_pos, int right_ref_pos, int left_diag, int right_diag, int left_col, int right_col, cost_t cost, std::vector<split_alignment_t>* target) const;
	/** Prepare an insertion-like instance of split_alignment_t and push it to a vector. */
	void push_insertion(const anchor_alignment_t& left_anchor, const anchor_alignment_t& right_anchor, int ref_breakpoint, int left_query_pos, int right_query_pos, int left_diag_idx, int right_diag_idx, int left_col_idx, int right_col_idx, cost_t cost, std::vector<split_alignment_t>* target) const;
	/** Prepare an interchromosomal rearrangement instance of split_alignment_t and push it to a vector. */
	void push_interchromosomal(const anchor_alignment_t& anchor1, const anchor_alignment_t& anchor2,  int breakpoint1_ref_pos, int breakpoint2_ref_pos, int diag1, int diag2, int col1, int col2, cost_t cost, std::vector<interchrom_split_alignment_t>* target) const;
public:
	SplitAligner(
		int cost_threshold = 115,
		gap_costs_t gap_costs = LINEAR,
		int indel_cost = 30,
		int gap_extend_cost = -1,
		int split_cost = 25,
		int interchrom_split_cost = 60,
		int min_anchor_length = 8,
		int secondary_aln_phred_diff = 29,
		int phred_offset = 33,
		bool simple_cigars = false,
		bool rightmost = false,
		bool allow_soft_clipping = false,
		int soft_clip_open_cost = 35,
		int soft_clip_extend_cost = 3
	);
	virtual ~SplitAligner();

	std::auto_ptr<std::vector<split_alignment_t> > combineAnchors(const anchor_alignment_t& left_anchor, const anchor_alignment_t& right_anchor) const;

	/** Tries to combine anchors from different chromosomes. Only returns alignments if not too much worse than the best known score given
	 *  as parameter. The maximal allowed difference is specified at construction time.
	 */
	std::auto_ptr<std::vector<interchrom_split_alignment_t> > combineAnchorsInterchromosomal(const anchor_alignment_t& left_anchor, const anchor_alignment_t& right_anchor, int best_known_score) const;

	/** Returns true if searching for an interchromosomal split alignment could yield a relevant alignment (given the score of the best known alignment). */
	bool doInterchromosomalSearch(int best_known_score) const;

	/** Tests whether an anchor alignment covers the whole read and, if so, returns a record
	 *  containing phred score and cigar string. */
	std::auto_ptr<whole_alignment_t> toWholeAlignment(const anchor_alignment_t& anchor) const;

	/** Computes a banded alignment between query and reference (starting from given position).
	 *  A variable length prefix of the query is used. The match is extended until the query is used completely or
	 *  the cost exceeds the cost given at construction time. Two such alignment anchors can later be combined into
	 *  one split read alignment.
	 *  If type is RIGHT, the query is matched from right to left against the reference, starting from ref_start.
	 */
	template <typename RefType, typename QueryType>
	std::auto_ptr<anchor_alignment_t> computeAnchorAlignment(const RefType& ref, const QueryType& query, size_t ref_start, anchor_type_t type) const {
		switch (gap_costs) {
			case LINEAR: 
				return computeAnchorAlignmentLinear(ref, query, ref_start, type);
			case AFFINE:
				return computeAnchorAlignmentAffine(ref, query, ref_start, type);
			default:
				assert(false);
		}
	}

	/** Print DP table of anchor alignment to given output stream. */
	template <typename RefType, typename QueryType, typename TableEntryType>
	void printTable(const anchor_alignment_t& anchor, const std::vector<std::vector<TableEntryType> >& table, std::ostream& out, const RefType& ref, const QueryType& query) const {
		int rows = anchor.column_count + max_indels;
		if (anchor.type == LEFT_ANCHOR) {
			rows = std::min(rows, (int)ref.size()-(int)anchor.ref_start_pos + 1);
		} else {
			rows = std::min(rows, (int)anchor.ref_start_pos + 2);
		}
		out << "                    - ";
		for (int j=1; j<(int)anchor.column_count; ++j) {
			out <<  std::setw(4) << j << " ";
		}
		out << "     " << std::endl;
		out << "                    - ";
		for (int j=1; j<(int)anchor.column_count; ++j) {
			if (anchor.type == LEFT_ANCHOR) {
				out <<  std::setw(4) << (query.qualityChar(j-1) - phred_offset) << " ";
			} else {
				out <<  std::setw(4) << (query.qualityChar(query.size()-j) - phred_offset) << " ";
			}
		}
		out << "     " << std::endl;
		out << "                    - ";
		for (int j=1; j<(int)anchor.column_count; ++j) {
			if (anchor.type == LEFT_ANCHOR) {
				out << "   " << query[j-1] << " ";
			} else {
				out << "   " << query[query.size()-j] << " ";
			}
		}
		out << " MIN " << std::endl;
		for (int i=0; i<rows; ++i) {
			out << std::setw(4) << i << " ";
			if (i == 0) {
				out << "          - ";
			} else {
				if (anchor.type == LEFT_ANCHOR) {
					out << std::setw(9) << (anchor.ref_start_pos + i - 1) << " ";
					out << ref[anchor.ref_start_pos + i - 1] << " ";
				} else {
					out << std::setw(9) << (anchor.ref_start_pos - i + 1) << " ";
					out << ref[anchor.ref_start_pos - i + 1] << " ";
				}
			}
			for (int j=0; j<(int)anchor.column_count; ++j) {
				int k = i - j + max_indels;
				if ((k>=0) && (k<(int)table.size())) {
					out << std::setw(4) << table[k][j] << " ";
				} else {
					out << "   - ";
				}
			}
			out << std::setw(4) << anchor.row_minima[i] << " " << std::endl;
		}
		out << " MIN         - - ";
		for (int j=0; j<(int)anchor.column_count; ++j) {
			out <<  std::setw(4) << anchor.column_minima[j] << " ";
		}
		out << "   - " << std::endl;
		out << std::endl;
	}

private:
	/** Anchor alignment core routine for linear gap costs.*/
	template <typename RefType, typename QueryType>
	std::auto_ptr<anchor_alignment_t> computeAnchorAlignmentLinear(const RefType& ref, const QueryType& query, size_t ref_start, anchor_type_t type) const {
		assert(ref_start < ref.size());
		assert(query.size() > 0);
		std::auto_ptr<anchor_alignment_t> result(new anchor_alignment_t(query.size(), ref_start, max_indels, type, cost_threshold, false));
		std::vector<std::vector<cost_t> >& d = result->diagonals;
		// Initialize first column
		for (int i=0; i<=max_indels; ++i) {
			int k = i + max_indels;
			d[k][0] = i * indel_cost;
			if (i>0) result->backtrace[k][0] = ABOVE;
		}
		for (int col = 1; col<=(int)query.size(); ++col) {
			cost_t minimum = std::numeric_limits<cost_t>::max();
			int min_row = std::max(0,col-max_indels);
			int max_row;
			if (type == RIGHT_ANCHOR) {
				max_row = std::min(col+max_indels, ((int)ref_start)+1);
			} else {
				max_row = std::min(col+max_indels,((int)ref.size())-((int)ref_start));
			}
			for (int row=min_row; row<=max_row; ++row) {
				// std::cout << "cell: " << row << ", " << col << std::endl;
				cost_t value = std::numeric_limits<cost_t>::max();
				backtrace_t backtrace = NONE;
				int diag = row - col + max_indels;
				// use different order depending on whether we have a left or right anchor
				// in order to report leftmost/rightmost indels
				if (type == LEFT_ANCHOR) { 
					if (rightmost) {
						// cell above
						if (diag>0) {
							cost_t c = d[diag-1][col] + indel_cost;
							if (c < value) { value = c; backtrace = ABOVE; }
						}
						// cell left
						if (diag+1<(int)d.size()) {
							cost_t c = d[diag+1][col-1] + indel_cost;
							if (c < value) { value = c; backtrace = LEFT; }
						}
						// cell up left
						if ((row>0) && (col>0)) {
							cost_t match = 0;
							// std::cout << "  comparing " << query[col-1] << " (" << (query.qualityChar(col-1) - phred_offset) << ") to " << ref[ref_start+row-1];
							if (query[col-1] != ref[ref_start+row-1]) {
								match = query.qualityChar(col-1) - phred_offset;
							}
							// std::cout << " --> match=" << match << std::endl;
							cost_t c = d[diag][col-1] + match;
							if (c < value) { value = c; backtrace = ABOVE_LEFT; }
						}
					} else {
						// cell up left
						if ((row>0) && (col>0)) {
							cost_t match = 0;
							// std::cout << "  comparing " << query[col-1] << " (" << (query.qualityChar(col-1) - phred_offset) << ") to " << ref[ref_start+row-1];
							if (query[col-1] != ref[ref_start+row-1]) {
								match = query.qualityChar(col-1) - phred_offset;
							}
							// std::cout << " --> match=" << match << std::endl;
							cost_t c = d[diag][col-1] + match;
							if (c < value) { value = c; backtrace = ABOVE_LEFT; }
						}
						// cell above
						if (diag>0) {
							cost_t c = d[diag-1][col] + indel_cost;
							if (c < value) { value = c; backtrace = ABOVE; }
						}
						// cell left
						if (diag+1<(int)d.size()) {
							cost_t c = d[diag+1][col-1] + indel_cost;
							if (c < value) { value = c; backtrace = LEFT; }
						}
					}
				} else {
					if (rightmost) {
						// cell up left
						if ((row>0) && (col>0)) {
							cost_t match = 0;
							if (query[query.size()-col] != ref[ref_start-row+1]) {
								match = query.qualityChar(query.size()-col) - phred_offset;
							}
							cost_t c = d[diag][col-1] + match;
							if (c < value) { value = c; backtrace = ABOVE_LEFT; }
						}
						// cell above
						if (diag>0) {
							cost_t c = d[diag-1][col] + indel_cost;
							if (c < value) { value = c; backtrace = ABOVE; }
						}
						// cell left
						if (diag+1<(int)d.size()) {
							cost_t c = d[diag+1][col-1] + indel_cost;
							if (c < value) { value = c; backtrace = LEFT; }
						}
					} else {
						// cell above
						if (diag>0) {
							cost_t c = d[diag-1][col] + indel_cost;
							if (c < value) { value = c; backtrace = ABOVE; }
						}
						// cell left
						if (diag+1<(int)d.size()) {
							cost_t c = d[diag+1][col-1] + indel_cost;
							if (c < value) { value = c; backtrace = LEFT; }
						}
						// cell up left
						if ((row>0) && (col>0)) {
							cost_t match = 0;
							if (query[query.size()-col] != ref[ref_start-row+1]) {
								match = query.qualityChar(query.size()-col) - phred_offset;
							}
							cost_t c = d[diag][col-1] + match;
							if (c < value) { value = c; backtrace = ABOVE_LEFT; }
						}
					}
				}
				d[diag][col] = value;
				result->backtrace[diag][col] = backtrace;
				result->row_minima[row] = std::min(result->row_minima[row], value);
				minimum = std::min(minimum, value);
			}
			// Test whether column contained any entries...
			if (minimum < std::numeric_limits<cost_t>::max()) {
				result->column_minima[col] = minimum;
				result->column_count = col + 1;
				result->row_count = max_row;
			} else {
				// ... if not, break (end of reference sequence reached)...
				break;
			}
			if (minimum > cost_threshold) {
				break;
			}
		}
		// printTable(*result, result->diagonals, std::cerr, ref, query);
		return result;
	}

	/** Anchor alignment core routine for affine gap costs.*/
	template <typename RefType, typename QueryType>
	std::auto_ptr<anchor_alignment_t> computeAnchorAlignmentAffine(const RefType& ref, const QueryType& query, size_t ref_start, anchor_type_t type) const {
		assert(ref_start < ref.size());
		assert(query.size() > 0);
		std::auto_ptr<anchor_alignment_t> result(new anchor_alignment_t(query.size(), ref_start, max_indels, type, cost_threshold, true));
		std::vector<std::vector<cost_t> >& d = result->diagonals;
		// Initialize first column
		for (int i=0; i<=max_indels; ++i) {
			int k = i + max_indels;
			d[k][0] = i * indel_cost;
			if (i>0) {
				result->backtrace[k][0] = ABOVE;
				result->indel_length_table[k][0] = i;
			}
		}
		// DP table (band) where each cell contains the optimal score of alignments ending on a deletion
		std::vector<std::vector<cost_t> > del_table(1+2*max_indels, std::vector<cost_t>(query.size()+1,cost_threshold+1));
		// corresponding tabel containing the length of this last deletion
		std::vector<std::vector<int> > del_length_table(1+2*max_indels, std::vector<cost_t>(query.size()+1,-1));
		// DP table (band) where each cell contains the optimal score of alignments ending on a insertion
		std::vector<std::vector<cost_t> > ins_table(1+2*max_indels, std::vector<cost_t>(query.size()+1,cost_threshold+1));
		// corresponding tabel containing the length of this last insertion
		std::vector<std::vector<int> > ins_length_table(1+2*max_indels, std::vector<cost_t>(query.size()+1,-1));
		for (int col = 1; col<=(int)query.size(); ++col) {
			cost_t minimum = cost_threshold + 1;
			int min_row = std::max(0,col-max_indels);
			int max_row;
			if (type == RIGHT_ANCHOR) {
				max_row = std::min(col+max_indels, ((int)ref_start)+1);
			} else {
				max_row = std::min(col+max_indels,((int)ref.size())-((int)ref_start));
			}
			for (int row=min_row; row<=max_row; ++row) {
				// std::cout << "cell: " << row << ", " << col << std::endl;
				cost_t value = std::numeric_limits<cost_t>::max();
				backtrace_t backtrace = NONE;
				int diag = row - col + max_indels;
				// compute entries for del_table and del_length_table
				cost_t del_cost = cost_threshold + 1;
				if (diag>0) {
					cost_t new_gap_cost = d[diag-1][col] + indel_cost + gap_extend_cost;
					cost_t extended_gap_cost = del_table[diag-1][col] + gap_extend_cost;
					if (new_gap_cost < extended_gap_cost) {
						del_cost = new_gap_cost;
						del_length_table[diag][col] = 1;
					} else {
						del_cost = extended_gap_cost;
						del_length_table[diag][col] = del_length_table[diag-1][col] + 1;
					}
					del_table[diag][col] = del_cost;
				}
				// compute entries for ins_table and ins_length_table
				cost_t ins_cost = cost_threshold + 1;
				if (diag+1<(int)d.size()) {
					cost_t new_gap_cost = d[diag+1][col-1] + indel_cost + gap_extend_cost;
					cost_t extended_gap_cost = ins_table[diag+1][col-1] + gap_extend_cost;
					if (new_gap_cost < extended_gap_cost) {
						ins_cost = new_gap_cost;
						ins_length_table[diag][col] = 1;
					} else {
						ins_cost = extended_gap_cost;
						ins_length_table[diag][col] = ins_length_table[diag+1][col-1] + 1;
					}
					ins_table[diag][col] = ins_cost;
				}
				// use different order depending on whether we have a left or right anchor
				// in order to report leftmost/rightmost indels
				if (type == LEFT_ANCHOR) {
					// compute cost that would incur from a match
					cost_t match_cost = cost_threshold + 1;
					if ((row>0) && (col>0)) {
						cost_t match = 0;
						// std::cout << "  comparing " << query[col-1] << " (" << (query.qualityChar(col-1) - phred_offset) << ") to " << ref[ref_start+row-1];
						if (query[col-1] != ref[ref_start+row-1]) {
							match = query.qualityChar(col-1) - phred_offset;
						}
						// std::cout << " --> match=" << match << std::endl;
						match_cost = d[diag][col-1] + match;
					}
					// update main DP table d
					if (rightmost) {
						// cell above
						if (del_cost < value) { value = del_cost; backtrace = ABOVE; }
						// cell left
						if (ins_cost < value) { value = ins_cost; backtrace = LEFT; }
						// cell up left
						if (match_cost < value) { value = match_cost; backtrace = ABOVE_LEFT; }
					} else {
						// cell up left
						if (match_cost < value) { value = match_cost; backtrace = ABOVE_LEFT; }
						// cell above
						if (del_cost < value) { value = del_cost; backtrace = ABOVE; }
						// cell left
						if (ins_cost < value) { value = ins_cost; backtrace = LEFT; }
					}
				} else {
					// compute cost that would incur from a match
					cost_t match_cost = cost_threshold + 1;
					if ((row>0) && (col>0)) {
						cost_t match = 0;
						if (query[query.size()-col] != ref[ref_start-row+1]) {
							match = query.qualityChar(query.size()-col) - phred_offset;
						}
						match_cost = d[diag][col-1] + match;
					}
					if (rightmost) {
						// cell up left
						if (match_cost < value) { value = match_cost; backtrace = ABOVE_LEFT; }
						// cell above
						if (del_cost < value) { value = del_cost; backtrace = ABOVE; }
						// cell left
						if (ins_cost < value) { value = ins_cost; backtrace = LEFT; }
					} else {
						// cell above
						if (del_cost < value) { value = del_cost; backtrace = ABOVE; }
						// cell left
						if (ins_cost < value) { value = ins_cost; backtrace = LEFT; }
						// cell up left
						if (match_cost < value) { value = match_cost; backtrace = ABOVE_LEFT; }
					}
				}
				d[diag][col] = value;
				result->backtrace[diag][col] = backtrace;
				result->row_minima[row] = std::min(result->row_minima[row], value);
				minimum = std::min(minimum, value);
				switch (backtrace) {
					case LEFT: result->indel_length_table[diag][col] = ins_length_table[diag][col]; break;
					case ABOVE: result->indel_length_table[diag][col] = del_length_table[diag][col]; break;
					case ABOVE_LEFT: break;
					default: assert(false);
				}
			}
			result->column_minima[col] = minimum;
			result->column_count = col + 1;
			result->row_count = max_row;
			if (minimum > cost_threshold) {
				break;
			}
		}
		// printTable(*result, result->diagonals, std::cerr, ref, query);
		return result;
	}
};

#endif /* SPLITALIGNER_H_ */
